library(phyloseq)
library(vegan)
library(dplyr)
library(reshape2)
library(here)
library(ggplot2)
library(DESeq2)
library(ggsci)
library(Biostrings)
library(seqinr)
library(janitor)
library(cowplot)
library(forcats)
library(tidyr)
library(DECIPHER)
library(phangorn)
library(ggrepel)
npg <- pal_d3("category10")(10)
theme_char <- function(base_size = 11, base_family = ""){
  theme_bw() %+replace%
    theme(axis.text = element_text(color = "Black"))
}
theme_set(theme_char())
set.seed(12)
ps <- readRDS(here("marcos_flea", "ps_trial.rds"))
sample_data <- read.csv(here("marcos_flea", "Sample_Info.csv"))
  rownames(sample_data) <- sample_data$ID
sample_data(ps) <- sample_data
psprop <- ps %>%
  transform_sample_counts(function(otu) otu/sum(otu))
multi <- genefilter_sample(ps, filterfun_sample(function(x) x > 0), A = 2)
multips <- prune_taxa(multi, ps)
  multips.glom <- tax_glom(psprop, taxrank = rank_names(ps)[6])
tax_table <- as.data.frame(tax_table(ps))
  tax_table$OTU <- rownames(tax_table)
seq <- refseq(ps)
ps
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 483 taxa and 10 samples ]
## sample_data() Sample Data:       [ 10 samples by 22 sample variables ]
## tax_table()   Taxonomy Table:    [ 483 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 483 reference sequences ]
ps.glom <- tax_glom(ps, taxrank = rank_names(ps)[6])
ps.melt <-  psmelt(ps)
psp.glom <- tax_glom(psprop, taxrank = rank_names(ps)[6])
  pspg.melt <- psmelt(psp.glom)
psp.melt <- psprop %>% psmelt
otu_table <- as.data.frame(otu_table(ps))
  otu_sums <- as.data.frame(colSums(otu_table))
  colnames(otu_sums) <- "Abundance"
  otu_sums$OTU <- rownames(otu_sums)
tax_table <- left_join(tax_table, otu_sums, by = "OTU")
wolb.melt <- psp.melt %>% filter(Genus == "Wolbachia")
  wolbotu <- unique(wolb.melt$OTU)
  wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV1", "ASV1", "Other")
  wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV2", "ASV2", wolb.melt$Big3)
  wolb.melt$Big3 <- ifelse(wolb.melt$OTU == "ASV3", "ASV3", wolb.melt$Big3)
wolb_otu_tab <- otu_table %>% select(wolbotu)
ps.wolb <- prune_taxa(wolbotu, ps)
bart.melt <- ps.melt %>% filter(Genus == "Bartonella")
  bartotu <- unique(bart.melt$OTU)
allfam <- names(sort(taxa_sums(psp.glom), decreasing=TRUE))
top20 <- names(sort(taxa_sums(psp.glom), decreasing=TRUE))[1:20]
ps.top20 <- transform_sample_counts(psp.glom, function(OTU) OTU/sum(OTU))
ps.top20 <- prune_taxa(top20, ps.top20)
plot_bar(ps.top20, x="ID", fill="Family")+ facet_wrap(~Time, scales="free_x")

ggsave(here("marcos_flea", "Image", "top20.png"), plot = last_plot(), width = 9, height = 6)
tax_table %>% dplyr::filter(Genus == "Wolbachia") %>% dplyr::mutate(sum = sum(Abundance))
##     Kingdom         Phylum               Class         Order          Family
## 1  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 2  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 3  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 4  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 5  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 6  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 7  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 8  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 9  Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 10 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 11 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 12 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 13 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 14 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 15 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 16 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 17 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 18 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 19 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 20 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 21 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 22 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 23 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 24 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 25 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 26 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 27 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 28 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 29 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 30 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 31 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 32 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 33 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
## 34 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales Anaplasmataceae
##        Genus Species     OTU Abundance     sum
## 1  Wolbachia    <NA>    ASV1   1336358 3159186
## 2  Wolbachia    <NA>    ASV2   1331247 3159186
## 3  Wolbachia    <NA>    ASV3    392785 3159186
## 4  Wolbachia    <NA>    ASV9     30562 3159186
## 5  Wolbachia    <NA>   ASV10     30379 3159186
## 6  Wolbachia    <NA>   ASV22     17357 3159186
## 7  Wolbachia    <NA>   ASV41      9373 3159186
## 8  Wolbachia    <NA>  ASV131      2705 3159186
## 9  Wolbachia    <NA>  ASV134      2673 3159186
## 10 Wolbachia    <NA>  ASV232      1324 3159186
## 11 Wolbachia    <NA>  ASV233      1319 3159186
## 12 Wolbachia    <NA>  ASV390       413 3159186
## 13 Wolbachia    <NA>  ASV450       298 3159186
## 14 Wolbachia    <NA>  ASV468       280 3159186
## 15 Wolbachia    <NA>  ASV502       246 3159186
## 16 Wolbachia    <NA>  ASV538       212 3159186
## 17 Wolbachia    <NA>  ASV564       191 3159186
## 18 Wolbachia    <NA>  ASV596       167 3159186
## 19 Wolbachia    <NA>  ASV599       164 3159186
## 20 Wolbachia    <NA>  ASV609       159 3159186
## 21 Wolbachia    <NA>  ASV617       157 3159186
## 22 Wolbachia    <NA>  ASV679       128 3159186
## 23 Wolbachia    <NA>  ASV714       113 3159186
## 24 Wolbachia    <NA>  ASV789        91 3159186
## 25 Wolbachia    <NA>  ASV844        71 3159186
## 26 Wolbachia    <NA>  ASV850        69 3159186
## 27 Wolbachia    <NA>  ASV864        62 3159186
## 28 Wolbachia    <NA>  ASV951        46 3159186
## 29 Wolbachia    <NA>  ASV960        45 3159186
## 30 Wolbachia    <NA>  ASV962        45 3159186
## 31 Wolbachia    <NA> ASV1011        38 3159186
## 32 Wolbachia    <NA> ASV1013        38 3159186
## 33 Wolbachia    <NA> ASV1021        37 3159186
## 34 Wolbachia    <NA> ASV1039        34 3159186
3159186/sum(otu_table)
## [1] 0.9666248
fam_replace <- tax_table
fam_replace$Family <- replace_na(fam_replace$Family, "Unassigned") 
fam_freq <- fam_replace %>% 
# mutate(Family = fct_infreq(Family)) %>% 
  ggplot(aes(x = Family))+
  geom_bar()+
  coord_flip()+
  labs(y = "Number of ASVs")
fam_abund <- fam_replace %>% group_by(Family) %>% mutate(fam_freq = sum(Abundance)) %>% 
  distinct(Family, .keep_all = TRUE) %>%  
  ggplot(aes(x = Family, y = fam_freq))+
  geom_col()+
  coord_flip()+
  labs(y = "Abundance", x = "")+theme(axis.text.y = element_blank())
plot_grid(fam_freq, fam_abund, rel_widths = c(1, 0.5))

ggsave(here("marcos_flea", "Image", "Fig2.eps"), width = 8, height = 12.5)
gen_glom <- fam_replace %>% group_by(Genus) %>% mutate(fam_freq = sum(Abundance)) %>% 
  distinct(Family, .keep_all = TRUE)
fam_glom <- fam_replace %>% group_by(Family) %>% mutate(fam_freq = sum(Abundance)) %>%   mutate(fam_asvs = n()) %>% distinct(Family, .keep_all = TRUE)


tax_glom <- as.data.frame(tax_table(ps.glom))
  tax_glom$OTU <- rownames(tax_glom)
prop_table <- as.data.frame(otu_table(ps.glom))
  prop_sums <- as.data.frame(colSums(prop_table))
  colnames(prop_sums) <- "Abundance"
  prop_sums$OTU <- rownames(prop_sums)
otu_prop <- prop_sums %>% mutate(Proportion = (round(Abundance/sum(Abundance)*100, 4)))
taxa_prop <- left_join(tax_glom, otu_prop, by = "OTU")
print_genus <- data.frame(Family = taxa_prop$Family,
                          Genus = taxa_prop$Genus,
                          Proportion = taxa_prop$Proportion)
write.csv(print_genus, here("marcos_flea", "genus_prop.csv"))

Diversity

invsimpson <- as.data.frame(vegan::diversity(otu_table, index = "invsimpson"))
  colnames(invsimpson) <- "invsimpson"
  invsimpson$ID <- rownames(invsimpson)
shannon <- as.data.frame(vegan::diversity(otu_table, index = "shannon"))
  colnames(shannon) <- "shannon"
  shannon$ID <- rownames(shannon)
spprich <- as.data.frame(specnumber(otu_table)) #calc species richness
  colnames(spprich) <- "spprich"
  spprich$ID <- rownames(spprich)
even <- as.data.frame(shannon[,1]/log(spprich[,1])) #calc Pielou's evenness
  colnames(even) <- "even"
  even$ID <- rownames(spprich)
diversity <- left_join(invsimpson, shannon, by = "ID")
diversity <- left_join(diversity, spprich, by = "ID")
diversity <- left_join(diversity, even, by = "ID")
diversity <- left_join(diversity, sample_data, by = "ID")

diversity$Cat_Bart <-  ifelse(diversity$Cat_Bart == "Bartonella Naive", "Bartonella Uninfected", diversity$Cat_Bart)
diversity$group <- ifelse(diversity$ID == "50UF", "50UF", diversity$Cat)
diversity$group <- ifelse(diversity$ID == "47UF", "47UF", diversity$group)
diversity$Cat_Bart <- ifelse(diversity$ID == "47UF", "Unfed Flea", diversity$Cat_Bart)
diversity$Cat_Bart <- ifelse(diversity$ID == "50UF", "Unfed Flea", diversity$Cat_Bart)
shannon.gg <- diversity %>% 
  ggplot(aes(x = as.factor(Time), y = shannon, color = Cat_Bart, group = group))+ 
  geom_point(alpha = 0.8)+geom_line()+
  labs(x = "Time (days)", y = "Shannon Index")+
  theme(axis.text.x = element_text(), legend.position = "none")+
  scale_color_manual(values = c("red", "black", "grey50"))
invsimpgg <- diversity %>% 
  ggplot(aes(x = as.factor(Time), y = invsimpson, color = Cat_Bart, group = group))+ 
  geom_point(alpha = 0.8)+geom_line()+
  labs(x = "Time (days)", y = "Inverse Simpson Index")+
  theme(axis.text.x = element_text(), legend.position = "none")+
  scale_color_manual(values = c("red", "black", "grey50"))
spprichgg <- diversity %>% 
  ggplot(aes(x = as.factor(Time), y = spprich, color = Cat_Bart, group = group))+ 
  geom_point(alpha = 0.8)+geom_line()+
  labs(x = "Time (days)", y = "Species Richness")+
  theme(axis.text.x = element_text(), legend.position = "none")+
  scale_color_manual(values = c("red", "black", "grey50"))
evennessgg <- diversity %>% 
  ggplot(aes(x = as.factor(Time), y =even, color = Cat_Bart, group = group))+ 
  geom_point(alpha = 0.8)+geom_line()+
  labs(x = "Time (days)", y = "Pielou's Evenness")+
  theme(axis.text.x = element_text(), legend.position = "none")+
  scale_color_manual(values = c("red", "black", "grey50"))

bart.type <- c(
  expression(paste(italic("Bartonella"), " Infected")),
  expression(paste(italic("Bartonella"), " Uninfected")),
  "Unfed Flea"
)

fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))

legend_plot <- diversity %>% 
  ggplot(aes(x = as.factor(Time), y =even, color = Cat_Bart, group = group))+
  geom_point()+geom_line()+
  labs(color = fill.tit)+
  scale_color_manual(values = c("red", "black", "grey50"), labels = bart.type)+
  theme(legend.text.align = 0)
legend <- get_legend(legend_plot)

ugh <- plot_grid(spprichgg, evennessgg, shannon.gg, legend, nrow = 1, labels = c("A", "B", "C"))
save_plot(plot = ugh, here("marcos_flea", "Image", "Fig3.jpeg"), base_width = 8, base_height = 3.5)
ugh

bhshan <- diversity %>% filter(Cat_Bart == "Bartonella Infected")
noshan <- diversity %>% filter(Cat_Bart == "Bartonella Uninfected") 
t.test(bhshan$shannon, noshan$shannon, paired = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  bhshan$shannon and noshan$shannon
## t = 1.403, df = 3.6663, p-value = 0.2394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3082353  0.8941888
## sample estimates:
## mean of x mean of y 
##  1.543527  1.250551
bhshan <- diversity %>% filter(Cat_Bart == "Bartonella Infected") %>% filter(Time == "1")
noshan <- diversity %>% filter(Cat_Bart == "Bartonella Uninfected")%>% filter(Time == "1")
t.test(bhshan$shannon, noshan$shannon)
## 
##  Welch Two Sample t-test
## 
## data:  bhshan$shannon and noshan$shannon
## t = 3.0834, df = 1.9442, p-value = 0.0943
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.2370271  1.3280859
## sample estimates:
## mean of x mean of y 
##  1.872497  1.326967
d1shan <- diversity %>% filter(Time == "1")
d9shan <- diversity %>% filter(Time == "9")
t.test(d1shan$shannon, d9shan$shannon, paired = T)
## 
##  Paired t-test
## 
## data:  d1shan$shannon and d9shan$shannon
## t = 2.5038, df = 3, p-value = 0.08742
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.1098851  0.9206561
## sample estimates:
## mean difference 
##       0.4053855
bray <- vegdist(otu_table, "bray")
bray <- melt(as.matrix(bray), varnames = c("row", "col"))
ggplot(bray, aes(x = row, y = col, fill = value))+
  geom_tile()

NMDS = metaMDS(otu_table, k = 2)
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.04473556 
## Run 1 stress 0.09181492 
## Run 2 stress 0.1015291 
## Run 3 stress 0.0506377 
## Run 4 stress 0.09411252 
## Run 5 stress 0.05063741 
## Run 6 stress 0.05063741 
## Run 7 stress 0.04473556 
## ... Procrustes: rmse 5.726807e-06  max resid 1.358413e-05 
## ... Similar to previous best
## Run 8 stress 0.0506373 
## Run 9 stress 0.04473556 
## ... Procrustes: rmse 3.930524e-06  max resid 5.979458e-06 
## ... Similar to previous best
## Run 10 stress 0.04473556 
## ... Procrustes: rmse 2.150821e-05  max resid 4.993466e-05 
## ... Similar to previous best
## Run 11 stress 0.04473556 
## ... Procrustes: rmse 5.359859e-06  max resid 9.428021e-06 
## ... Similar to previous best
## Run 12 stress 0.04473556 
## ... Procrustes: rmse 2.865025e-06  max resid 6.84322e-06 
## ... Similar to previous best
## Run 13 stress 0.04473556 
## ... Procrustes: rmse 1.628114e-05  max resid 3.832866e-05 
## ... Similar to previous best
## Run 14 stress 0.04473556 
## ... Procrustes: rmse 2.123653e-06  max resid 3.700492e-06 
## ... Similar to previous best
## Run 15 stress 0.04473556 
## ... Procrustes: rmse 2.346207e-05  max resid 5.489858e-05 
## ... Similar to previous best
## Run 16 stress 0.05063741 
## Run 17 stress 0.08960182 
## Run 18 stress 0.05063743 
## Run 19 stress 0.09411241 
## Run 20 stress 0.09181554 
## *** Best solution repeated 8 times
sample.scores <- as.data.frame(scores(NMDS)[1])
  sample.scores$ID <- rownames(sample.scores)
  sample.scores <- left_join(sample.scores, sample_data, by = "ID")
asv.scores <- as.data.frame(scores(NMDS)[2])
sample.scores$xnudge <- "0.05"
sample.scores$ynudge <- "0.05"
sample.scores$xnudge <- ifelse(sample.scores$ID == "3711-24hr", 0.10, sample.scores$xnudge)
sample.scores$Cat_Bart <- ifelse(sample.scores$Cat_Bart == "Bartonella Uninfected", "Bartonella Uninfected", sample.scores$Cat_Bart)
sample.scores$ID <- ifelse(sample.scores$ID == "50UF", "UF50", sample.scores$ID)
sample.scores$ID <- ifelse(sample.scores$ID == "47UF", "UF47", sample.scores$ID)
sample.scores$Cat_Bart <- ifelse(sample.scores$ID == "UF47", "Unfed Flea", sample.scores$Cat_Bart)
sample.scores$Cat_Bart <- ifelse(sample.scores$ID == "UF50", "Unfed Flea", sample.scores$Cat_Bart)

bart.type <- c(
  expression(paste(italic("Bartonella"), " Infected")),
  expression(paste(italic("Bartonella"), " Uninfected")),
  "Unfed Flea"
)

fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))

ggplot(sample.scores, aes(x = sites.NMDS1, y = sites.NMDS2, color = Cat_Bart))+
  geom_point()+
  theme_char()+
  geom_text(aes(label = ID), nudge_x = 0.23)+
  labs(color = fill.tit)+
  guides(color = guide_legend(override.aes=list(shape = 15, size = 5)))+
  scale_color_manual(values = c("red", "black", "grey45"), labels = bart.type)+
  theme(legend.text.align = 0)

ggsave(here("marcos_flea", "Image", "Fig4.eps"), plot = last_plot(), width = 7, height = 4.5)

Bartonella

bart.melt %>% 
ggplot(aes(x = ID, y = Abundance, fill = Cat_Bart))+
  geom_col(color = "black")+theme_char()+
  facet_grid(OTU~Cat_Bart, scales = "free_x")+
  theme(axis.text.x = element_text(angle = 45,hjust = 1), legend.position = "none")+
  labs(y = "Abundance")+
  geom_text(aes(label = Abundance), vjust = -0.30)

ggsave(here("marcos_flea", "Image", "bart.png"), plot = last_plot(), width = 9, height = 6)

Wolbachia

wolbseq <- seq[wolbotu]
writeXStringSet(wolbseq, here("marcos_flea", "sequences", "wolbseq.fasta"))
pspg.melt %>% filter(Genus == "Wolbachia") %>% 
  ggplot(aes(y = Abundance, x = ID, fill = as.factor(Cat_Bart)))+
  geom_col(color = "black")+theme_char()+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(y = "Wolbachia Proportion", color = "Cat Bartonella")+
  facet_grid(~Time, scales = "free_x")+
  scale_fill_manual(values = c("Pink", "Black"))

wolb.melt %>% mutate(Cat = factor(Cat, levels = c("Unfed", "3363", "3508", "3320", "3711"))) %>% 
  ggplot(aes(x = ID, y = Abundance, fill = Big3))+
  geom_col(color = "black")+theme_char()+
  facet_wrap(~Cat, scales = "free_x", nrow = 1)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  labs(y = "Proportion of Total Reads")

ggsave(here("marcos_flea", "Image", "wolbbig3.png"), plot = last_plot(), width = 9, height = 6)
for(i in 1:10){
graph <- wolb.melt %>% filter(OTU == wolbotu[i]) %>% filter(Abundance > 0) %>% 
  ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat, group = Cat))+
  geom_point()+geom_line()+
  facet_wrap(~Cat_Bart, scales = "free_x")+theme_char()+
  labs(x = "Time (days)", y = "Proportion", title = paste(wolbotu[i], "Proportion"))
print(graph)
}

wolb.interest <- c("ASV1", "ASV2", "ASV3")

wolb.melt$group <- ifelse(wolb.melt$Sample == "50UF", "50UF", wolb.melt$Cat)
wolb.melt$group <- ifelse(wolb.melt$Sample == "47UF", "47UF", wolb.melt$group)
wolb.melt$group <- ifelse(wolb.melt$Sample == "47UF", "47UF", wolb.melt$group)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Sample == "47UF", "Unfed Flea", wolb.melt$Cat_Bart)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Sample == "50UF", "Unfed Flea", wolb.melt$Cat_Bart)
wolb.melt$Cat_Bart <- ifelse(wolb.melt$Cat_Bart == "Bartonella Naive", "Bartonella Uninfected", wolb.melt$Cat_Bart)

bart.type <- c(
  expression(paste(italic("Bartonella"), " Infected")),
  expression(paste(italic("Bartonella"), " Uninfected")),
  "Unfed Flea"
)

fill.tit <- expression(paste("Cat ", italic("Bartonella"), " Status"))
wolb.melt %>% filter(OTU %in% wolb.interest) %>% mutate()  %>% 
  mutate(OTU = factor(OTU, levels = c(wolb.interest))) %>% 
  mutate(Cat_Bart = factor(Cat_Bart, levels = c("Bartonella Infected", "Bartonella Uninfected", "Unfed Flea"))) %>% 
  ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart, group = group))+
  geom_line()+geom_point()+
  facet_wrap(~OTU)+theme_char()+
  labs(x = "Time (days)", y = "Proportion", color = fill.tit)+
  scale_color_manual(values = c("red","black", "grey65"), labels = bart.type)+
  theme(legend.position = "right", legend.text.align = 0)

ggsave(here("marcos_flea", "Image", "Fig5.tiff"), plot = last_plot(), width = 7, height = 3.5)
asv.interest <- c("ASV1", "ASV2", "ASV3", "ASV9", "ASV10", "ASV352")

psp.melt$group <- ifelse(psp.melt$Sample == "50UF", "50UF", psp.melt$Cat)
psp.melt$group <- ifelse(psp.melt$Sample == "47UF", "47UF", psp.melt$group)
psp.melt$group <- ifelse(psp.melt$Sample == "47UF", "47UF", psp.melt$group)
psp.melt$Cat_Bart <- ifelse(psp.melt$Sample == "47UF", "Unfed Flea", psp.melt$Cat_Bart)
psp.melt$Cat_Bart <- ifelse(psp.melt$Sample == "50UF", "Unfed Flea", psp.melt$Cat_Bart)
psp.melt %>% filter(OTU %in% asv.interest) %>% mutate()  %>% 
  mutate(OTU = factor(OTU, levels = c(asv.interest))) %>% 
  mutate(Cat_Bart = factor(Cat_Bart, levels = c("Bartonella Naive", "Bartonella Infected", "Unfed Flea"))) %>% 
  ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart, group = group))+
  geom_line()+geom_point()+
  facet_wrap(~OTU, scales = "free_y")+theme_char()+
  labs(x = "Time (days)", y = "Proportion", color = "Host Bartonella Status")+
  scale_color_manual(values = c("black", "red", "grey69"))+
  theme(legend.position = "bottom")

hm <-prune_taxa(names(sort(taxa_sums(ps),TRUE)[1:300]), ps)
plot_heatmap(hm, sample.label="ID", na.value = "white")

DeSeq

Cat Bart w/o nonfed

multinounfed <- subset_samples(multips, Fed == "Yes")
deseq <- phyloseq_to_deseq2(multinounfed, ~ Cat_Bart)
## converting counts to integer mode
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]
deseq <- DESeq(deseq, test = "Wald", fitType = "parametric")
## estimating size factors
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]
## final dispersion estimates
## fitting model and testing
res <- results(deseq, cooksCutoff = FALSE)
alpha = 0.01
res = cbind(as(res, "data.frame"), as(tax_table(multips)[rownames(res), ], "matrix"))
wolb.dseq <- res %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.dseq
##  [1] baseMean       log2FoldChange lfcSE          stat           pvalue        
##  [6] padj           Kingdom        Phylum         Class          Order         
## [11] Family         Genus          Species       
## <0 rows> (or 0-length row.names)
sigtab = res[which(res$padj < alpha),]
head(sigtab)
##         baseMean log2FoldChange    lfcSE      stat       pvalue         padj
## ASV288 127.53884      -24.79779 3.098097 -8.004200 1.202456e-15 1.394849e-13
## ASV336  77.19327      -24.11462 3.383905 -7.126269 1.031263e-12 3.987551e-11
## ASV377  59.77222      -23.76011 3.384060 -7.021185 2.199948e-12 5.103880e-11
## ASV391  76.36428      -23.88336 3.383911 -7.057917 1.690174e-12 4.901504e-11
## ASV402  31.38034      -22.89174 3.384681 -6.763339 1.348472e-11 2.214814e-10
## ASV432  24.36797      -21.79153 3.385056 -6.437570 1.214016e-10 1.280235e-09
##         Kingdom           Phylum          Class            Order
## ASV288 Bacteria       Firmicutes     Clostridia   Lachnospirales
## ASV336 Bacteria       Firmicutes     Clostridia   Lachnospirales
## ASV377 Bacteria       Firmicutes     Clostridia   Lachnospirales
## ASV391 Bacteria       Firmicutes     Clostridia   Lachnospirales
## ASV402 Bacteria       Firmicutes     Clostridia   Lachnospirales
## ASV432 Bacteria Actinobacteriota Coriobacteriia Coriobacteriales
##                 Family                           Genus          Species
## ASV288 Lachnospiraceae                    Oribacterium asaccharolyticum
## ASV336 Lachnospiraceae    [Ruminococcus] torques group             <NA>
## ASV377 Lachnospiraceae                            <NA>             <NA>
## ASV391 Lachnospiraceae                         Blautia        stercoris
## ASV402 Lachnospiraceae [Eubacterium] fissicatena group             <NA>
## ASV432 Eggerthellaceae                     Eggerthella            lenta
# Phylum order
x = tapply(sigtab$log2FoldChange, sigtab$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtab$Phylum = factor(as.character(sigtab$Phylum), levels=names(x))
# Genus order
x = tapply(sigtab$log2FoldChange, sigtab$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtab$Genus = factor(as.character(sigtab$Genus), levels=names(x))
sigtab$OTU <- rownames(sigtab)
sigtab$Analysis <- c("Cat Bartonella Status")
ggplot(sigtab, aes(x=Genus, y=log2FoldChange, color=Phylum)) + 
  geom_point(size=2) + 
  theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5))+coord_flip()+
  labs(title = "Cat Bartonella Status")

Time w/0 unfed

deseq.time <- phyloseq_to_deseq2(multinounfed, ~ Time)
## converting counts to integer mode
##   the design formula contains one or more numeric variables with integer values,
##   specifying a model with increasing fold change for higher values.
##   did you mean for this to be a factor? if so, first convert
##   this variable to a factor using the factor() function
deseq.time <- DESeq(deseq.time, test = "Wald", fitType = "parametric")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res.time <- results(deseq.time, cooksCutoff = FALSE)
alpha = 0.05
res.time = cbind(as(res.time, "data.frame"), as(tax_table(multinounfed)[rownames(res.time), ], "matrix"))
wolb.time <- res.time %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.time
##  [1] baseMean       log2FoldChange lfcSE          stat           pvalue        
##  [6] padj           Kingdom        Phylum         Class          Order         
## [11] Family         Genus          Species       
## <0 rows> (or 0-length row.names)
sigtab.time = res.time[which(res.time$padj < alpha),]
sigtab.time$OTU <- rownames(sigtab.time)
head(sigtab.time)
##         baseMean log2FoldChange     lfcSE      stat       pvalue         padj
## ASV193 145.90173      -3.563588 0.3815665 -9.339363 9.691577e-21 1.124223e-18
## ASV279  68.48759      -3.428953 0.3789772 -9.047915 1.457247e-19 8.452035e-18
## ASV314  75.77528      -3.422823 0.3828236 -8.940993 3.856871e-19 1.223514e-17
## ASV336  77.19327       2.468604 0.4229935  5.836034 5.345796e-09 3.875702e-08
## ASV367  68.10692      -3.424755 0.3834652 -8.931072 4.219012e-19 1.223514e-17
## ASV377  59.77222       2.985442 0.4230134  7.057559 1.694531e-12 1.310438e-11
##         Kingdom     Phylum      Class           Order          Family
## ASV193 Bacteria Firmicutes Clostridia Oscillospirales Ruminococcaceae
## ASV279 Bacteria Firmicutes Clostridia  Lachnospirales Lachnospiraceae
## ASV314 Bacteria Firmicutes Clostridia  Lachnospirales Lachnospiraceae
## ASV336 Bacteria Firmicutes Clostridia  Lachnospirales Lachnospiraceae
## ASV367 Bacteria Firmicutes Clostridia  Lachnospirales Lachnospiraceae
## ASV377 Bacteria Firmicutes Clostridia  Lachnospirales Lachnospiraceae
##                                Genus Species    OTU
## ASV193                Incertae Sedis    <NA> ASV193
## ASV279             Lachnoclostridium    <NA> ASV279
## ASV314                          <NA>    <NA> ASV314
## ASV336  [Ruminococcus] torques group    <NA> ASV336
## ASV367 Lachnospiraceae NK4A136 group    <NA> ASV367
## ASV377                          <NA>    <NA> ASV377
# Phylum order
x.time = tapply(sigtab.time$log2FoldChange, sigtab.time$Phylum, function(x) max(x))
x.time = sort(x.time, TRUE)
sigtab.time$Phylum = factor(as.character(sigtab.time$Phylum), levels=names(x.time))
# Genus order
x.time = tapply(sigtab.time$log2FoldChange, sigtab.time$Genus, function(x) max(x.time))
x.time = sort(x.time, TRUE)
sigtab.time$Genus = factor(as.character(sigtab.time$Genus), levels=names(x.time))
sigtab.time$Analysis <- c("Time")
sigtab.time$OTU <- rownames(sigtab.time)
ggplot(sigtab.time, aes(x = 1, y = OTU, fill = log2FoldChange))+
  geom_tile()+
  theme(axis.text.y = element_text(size = 6), axis.text.x = element_blank())+
  facet_wrap(~Phylum, scales = "free_y")

ggplot(sigtab.time, aes(x=Genus, y=log2FoldChange, color=Phylum)) + 
  geom_point( alpha = 0.8) + 
  theme(axis.text.x = element_text(),
        legend.position = "bottom")+coord_flip()

Fed Status

deseq.fed <- phyloseq_to_deseq2(multips, ~ Fed)
## converting counts to integer mode
deseq.fed <- DESeq(deseq.fed, test = "Wald", fitType = "parametric")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## -- replacing outliers and refitting for 43 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)
## estimating dispersions
## fitting model and testing
res.fed <- results(deseq.fed, cooksCutoff = FALSE)
res.fed = cbind(as(res.fed, "data.frame"), as(tax_table(multips)[rownames(res.fed), ], "matrix"))
wolb.fed <- res.fed %>% filter(Genus == "Wolbachia") %>% filter(padj < 0.05)
wolb.fed 
##        baseMean log2FoldChange   lfcSE     stat       pvalue         padj
## ASV390 12.35088       18.66602 3.88657 4.802698 1.565417e-06 3.506535e-05
##         Kingdom         Phylum               Class         Order
## ASV390 Bacteria Proteobacteria Alphaproteobacteria Rickettsiales
##                 Family     Genus Species
## ASV390 Anaplasmataceae Wolbachia    <NA>
sigtab.fed = res.fed[which(res.fed$padj < alpha),]
alpha = 0.01
sigtab.fed$OTU <- rownames(sigtab.fed)
head(sigtab.fed)
##        baseMean log2FoldChange    lfcSE     stat       pvalue         padj
## ASV285 23.91761       19.57942 3.885561 5.039019 4.679249e-07 1.903797e-05
## ASV314 24.58057       19.61182 3.885533 5.047394 4.478771e-07 1.903797e-05
## ASV367 14.15733       18.85508 3.886304 4.851674 1.224237e-06 3.427863e-05
## ASV390 12.35088       18.66602 3.886570 4.802698 1.565417e-06 3.506535e-05
## ASV400 14.36153       17.76983 3.886279 4.572453 4.820462e-06 7.712740e-05
## ASV403 23.06325       19.51557 3.885603 5.022533 5.099456e-07 1.903797e-05
##         Kingdom          Phylum               Class            Order
## ASV285 Bacteria      Firmicutes          Clostridia  Oscillospirales
## ASV314 Bacteria      Firmicutes          Clostridia   Lachnospirales
## ASV367 Bacteria      Firmicutes          Clostridia   Lachnospirales
## ASV390 Bacteria  Proteobacteria Alphaproteobacteria    Rickettsiales
## ASV400 Bacteria      Firmicutes          Clostridia   Lachnospirales
## ASV403 Bacteria Acidobacteriota      Acidobacteriae Acidobacteriales
##                 Family                         Genus Species    OTU
## ASV285 Ruminococcaceae                Incertae Sedis    <NA> ASV285
## ASV314 Lachnospiraceae                          <NA>    <NA> ASV314
## ASV367 Lachnospiraceae Lachnospiraceae NK4A136 group    <NA> ASV367
## ASV390 Anaplasmataceae                     Wolbachia    <NA> ASV390
## ASV400 Lachnospiraceae                 GCA-900066575    <NA> ASV400
## ASV403            <NA>                          <NA>    <NA> ASV403
# Phylum order
x.fed = tapply(sigtab.fed$log2FoldChange, sigtab.fed$Phylum, function(x) max(x))
x.fed = sort(x.fed, TRUE)
sigtab.fed$Phylum = factor(as.character(sigtab.fed$Phylum), levels=names(x.fed))
# Genus order
x.fed = tapply(sigtab.fed$log2FoldChange, sigtab.fed$Genus, function(x) max(x.fed))
x.fed = sort(x.fed, TRUE)
sigtab.fed$Genus = factor(as.character(sigtab.fed$Genus), levels=names(x.fed))
sigtab.fed$Analysis <- c("Fed Status")

ggplot(sigtab.fed, aes(x = 1, y = OTU, fill = log2FoldChange))+
  geom_tile()+
  theme(axis.text.y = element_text(size = 6), axis.text.x = element_blank())+
  facet_wrap(~Phylum, scales = "free_y")

ggplot(sigtab.fed, aes(x=Genus, y=log2FoldChange, color=Phylum)) + 
  geom_point( alpha = 0.8) + 
  theme(axis.text.x = element_text(),
        legend.position = "right")+coord_flip()+
  labs(title = "Fed Status")

All Together

alldseq <- bind_rows(sigtab, sigtab.fed, sigtab.time)
alldseq <- alldseq[complete.cases(alldseq$Genus),]



ggplot(alldseq, aes(x=OTU, y=Analysis, fill = log2FoldChange)) + 
  geom_tile()+
  theme(axis.text.x = element_text(),
        legend.position = "right",
        axis.text.y = element_text(size = 10))+coord_flip()

ggplot(alldseq, aes(x=Genus, y=log2FoldChange, fill = Phylum)) + 
  geom_col(position = "dodge")+
  theme(axis.text.x = element_text(),
        legend.position = "bottom",
        axis.text.y = element_text(size = 10))+coord_flip()+
  facet_wrap(~Analysis)

alldseq %>% mutate(Analysis = factor(Analysis, levels = c("Fed Status", "Cat Bartonella Status", "Time"))) %>% 
  mutate(OTU = paste0(Genus, " (",OTU, ")")) %>% 
ggplot(aes(x=OTU, y=log2FoldChange)) + 
  geom_segment(aes(x = OTU, y = 0, xend = OTU, yend = log2FoldChange))+
  geom_point(size = 2, aes(color = Order))+
  theme(axis.text.x = element_text(),
        legend.position = "bottom",
        axis.text.y = element_text(size = 10))+
  coord_flip()+
  facet_wrap(~Analysis)

ggsave(here("marcos_flea", "Image", "alldeseq.png"), plot = last_plot(), width = 8, height = 5.5)
psp.melt %>% filter(OTU %in% alldseq$OTU) %>% filter(Abundance > 0) %>% 
  mutate(OTU = paste0(Genus, " (",OTU, ")")) %>% 
  ggplot(aes(x = as.factor(Time), y = Abundance, color = Cat_Bart))+
  geom_point()+
  facet_wrap(~OTU, nrow = 5)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "bottom", strip.text = element_text(size = 11))+
  labs(x = "Time (Days)", y = "Proportion of Reads")+
  scale_color_manual(values = c("Red", "Black", "grey"))

ggsave(here("marcos_flea", "Image", "individualdeseq.png"), plot = last_plot(), width = 13, height = 7)
psprop %>% psmelt %>% filter(Genus == "Blautia") %>%  filter(Abundance > 0) %>% 
  ggplot(aes(y = Abundance, x = as.factor(Time), fill = Cat_Bart))+
  geom_col(aes(color = Cat_Bart))+geom_point(color = "Black", alpha = 0.6)+
  theme(axis.text.x = element_text(angle = 90, hjust = 1), legend.position = "right")+
  labs(x = "Time (Days)", y = "Proportion of Reads", title = "Blautia Abundance")